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Positive deconvolution for superimposed 
extended source and point sources. 
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Abstract. The paper deals with the construction of images from visibilities acquired using aperture synthesis 
instruments: Fourier synthesis, deconvolution, and spectral interpolation/extrapolation. Its intended application 
is to specific situations in which the imaged object possesses two superimposed components: (i) an extended 
component together with (ii) a set of point sources. It is also specifically designed to the case of positive maps, and 
accounts for a known support. Its originality lies within joint estimation of the two components, coherently with 
data, properties of each component, positivity and possible support. We approach the subject as an inverse problem 
within a regularization framework: a regularized least-squares criterion is specifically proposed and the estimated 
maps are defined as its minimizer. We have investigated several options for the numerical minimization and we 
propose a new efficient algorithm based on augmented Lagrangian. Evaluation is carried out using simulated and 
real data (from radio interferometry) demonstrating the capability to accurately separate the two components. 

Key words. Techniques: image processing, interferometric, deconvolution, inverse problem, extended/point sources. 



1. Introduction 

Radio interferometers can be seen as instruments mea- 
suring a set of 2D-Fourier coefficients (visibilities) of the 
brightness distribution of a region in the sky. Visibilities 
are measured in the Fourier domain (the (u, v)-plane) by 
means of different baselines (projected distance between 
cross-correlated a ntennas) . Practically, t here are two prin- 
cipal deficiencies llThompson et al . 2001) in the visibilities 

1. the limited coverage of the (it, w)-plane, 

2. measurements errors (especially in millimeter range). 

Regarding point 1, three limitations are encountered. 

— Usually the central part of the aperture (up to the 
antenna diameter) is not observed. From this stand 
point, interferometers behave as high pass filters. 

— Information above the longest baseline is unavailable. 
In this sense the instruments behave as low pass filters. 

— The (u, u)-plane coverage is irregular, especially when 
there is a small number of antennas. This results in 
dirty beam (Fourier transform of visibility weights) 
with intricate structure and strong sidelobes. 

Thus, such instruments can be seen as band pass fil- 
ters with an intricate impulse response (dirty beam), and 
noisy output. As a consequence, the available data is rela- 
tively poor for imaging objects with various spatial struc- 



tures extended over the whole frequency domain. In or- 
der to compensate for these deficiencies, a large number 
of methods (from model fitting to non parametric de- 
co nvolution) has been continuously proposed (see review 
m llStarck et alJ l200211 and specialized for different types 
of maps. The present paper deals with a particular type of 
map consisting of the superimposition of two components. 

— Point Source (PS), or nearly black objects: essentially 
null-component, with a few strong point sources. 

— Extended Sources (ES): spatially extended, smooth 
components. 

The problem at hand is to build reliable and accurate es- 
timates of two distinct maps (one for PS, one for ES) from 
a unique given set of visibilities. The question arises e.g. 
for radio imaging of the solar corona at meter wavelength 
where very strong storms are superimposed over a more 
stable and large quiet Sun radio-emission (see Sect. 14.11) . 

Remark 1. From a statistical standpoint, PS/ES can be 
modelized as set of uncorrelated/correlated pixels, respec- 
tively. In the Fourier plane they are respectively charac- 
terized by an extension over the whole frequency domain 
(PS) and an extension reduced to the low frequencies do- 
main (ES). In particular, both of them have significant 
components in the low frequencies domain (see Fig. ^-a) . 
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Fig. 1. a: The solid line with squares (resp. dashed line 
with circles) shows spectral content for ES (resp. PS). 
Both of them have low frequencies components, b: The 
two lines show spectral contents of correlated components 
(ES), with different level of correlation, c: Elementary de- 
composition for wavelet transform. The solid line with 
squares (resp. dashed line with circles) shows low (resp. 
high) frequency content. 

1.1. General bibliographical analysis 

In order to compensate for the deficiencies in the avail- 
able data, additional information is (implicitly or explic- 
itly) accounted for. Practically, most existing methods are 
founded on specific expected properties of observed and 
reconstructed sources. The proposed analysis relies on un- 
derlying decompositions of unknown image. 

PS based methods - A first part of existing methods 
relies on PS prop erties. Into this category fall or iginal ver- 
sions of CLEAN <|HogbonJl974llFomalontll973|) . which it- 
eratively withdraw the PS contribution to th e dirty map . 
Early Maximum Entropy Methods (MEM) (jAbleslll97^ 
are also founded on the properties of PS: in a regular- 
ized context, they introduce separable penalization terms 
(without pixel interaction) and favor high-amplitude PS. 

ES based methods - Two main classes of methods have 
been proposed to account for the correlation of ES. 

— The correlation structure is introduced by a convolu- 
tion kernel. This is the case in MEM with an Intrinsic 
Correlati on Function (ICF ) llGullI Il989h and Pix on 
methods l|Dixon et al.lll996t IPuetter fc Yahillll999|) . 

— The other class of method relies on pixel interactive 
pena lty. The early versions inv olve quadratic penal- 
ties ijTikhonov fc Arseninl Il977h . Extensions to other 



pena l ties have also been widely developed JO'Sullivanl 
ll995llSnVder et al.lll992HMugnier et alJl2004l) . ~ 

Mixed ES+SP model - The case of an explicit model 
mixing ES and PS has also been addressed; however, 
literature in this case is poor. To our knowledge , two 
papers have been published: (|Magain et alJ ll99Sl) and 
llPirzkal et aTl boom. They introduced the decomposition 
of the search map as the sum of a PS map and an ES map. 
From a spectral standpoint, PS / ES are respectively char- 
acterized as shown in Fig. ^a (see also Remark . The 
present paper is founded on this approach (see Sect. I1.2|l . 

Multi-resolution / subband methods - Another class 
of method received a large attention, namely the multi- 
resolution and subband approaches. 



The a pproach proposed by ljWeirll992HBontekoe et all 
11994 introduces structure by means of different ICF. 
The unknown map is the sum of several ES, with dif- 
ferent level of correlation, i.e. several low frequency 
components. The underlying decomposition is shown 
in Fig.^b in the case of two components. 
We also have witnessed the dev elopment of multi- 
resolu tion extensions of CLEAN (Wa kker fc Schward 
1988) as well as more subtle appro aches based 
on wavelet decomposition and MEM llStarck et all 
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ll994l:IPantin fc Starckll996tlStarck et all200l|) . These 
methods are less specific and widely used for general 
deconvolution. They aim at reconstructing maps with 
different scales by splitting the Fourier plane into vari- 
ous zones. They basically rely on (recursive) decompo- 
sition in low and high frequencies as shown in Fig.^-c. 

1.2. PS plus ES: proposed developments 

As mentioned above, the present paper is devoted to the 
estimation of two distinct maps (one for ES and one for 
PS) from a unique set of visibilities. We the n naturally re- 
sort to the work of lMagain et alJ l|l998ft and lPirzkal et alJ 
(2000). In both cases the PS map is written in a para- 
metric manner founded on positions and amplitudes of 
peaks. Smoothness of the ES is in cluded by means o f 
Gaussian ICF and M EM penalty JPirzkal et all Eoooj) 
and Tikhonov penalty ijMagain et al.lfl998|l . Nevertheless. 
they both have sev eral limitations. On the one hand, 
l|Pirzkal et alJ 12000) relics on the knowledge of the po- 
sition of the PS which i s not available to us . On the other 
hand, the drawback of l|Magain et alJll99^) is twofold. 

1. It does not deconvolve with the total PSF. 

2. The optimized criterion is intricate w.r.t. the PS po- 
sitions so, it is not alway s possible to find th e global 
minimum of the criterion l|Magain et al.lll998l p. 474). 

On the contrary, our approach achieves a complete 
deconvolution. Moreover, our work introduces proper- 
ties so that an optimal solution is properly defined and 
practically attainable. In a unique coherent framework, 
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the proposed method simultaneously accounts for intri- 
cate dirty beam, noise, the existence of point sources 
superimposed onto a smooth component, positivity, and 
the possible knowledge of a support. The estimated 
maps are defined as the constrained minimizer of a 
penalized least-squares criterion specifically adapted to 
this situation. So that, the method assigns a coher- 
ent value to unmeasured Fourier coefficients. The basic 
ideas developed here hav e already been par tly presented 
with in spectral analysis Jciuciu et alJl2f)f)lJ) . spectrome- 
try (Mohammad-Dia fari et al.ll2002j) and satellite imag- 
ing ijsamson et al.ll2003|) . 

The paper is organized as follows. In Sect.El we define 
notations and state the problem in three classical forms: 
Fourier synthesis, spectral extrapolation / interpolation 
and deconvolution. All three cases concern rank-deficient 
linear inverse problems with additive noise. The proposed 
method is presented in Sect. [31 Section |3~T1 introduces the 
regularization principles used in the subsequent sections; 
Sections 13.21 and 13.31 respectively deal with PS and ES 
map; Section [3.41 is devoted to the main contribution: the 
reconstruction of two maps simultaneously, one consist- 
ing of PS, and the other of ES. Simulation and real data 
computations are presented throughout Sect. 01 From a 
numerical optimization viewpoint, the proposed method 
reduces to a constrained quadratic programming problem 
and various options have been studied and compared. The 
proposed algorithm founded on augmented Lagrangian 
principle is presented in Sect. In Sect. El we set out 
conclusions and perspectives. 

2. Problem statement and least squares solution 

The usual model 1 for the instrument writes as a weighted 
truncated noisy Fourier transform (discrete and regular): 

y = WTFx + b, (1) 

where x £ M> N is the unknown map and y and b £ <C M 
are the Fourier coefficient and noise (N unkown param- 
eters for M measurements). F is the N x N normalized 
FFT matrix and T is a 0/1 -binary truncation (or sam- 
pling) M x N matrix (T discards frequencies outside the 
(u, w)-plane coverage). W is a M x M diagonal matrix ac- 
counting for visibility weights. For the sake of simplicity 
and in accordance with real data processed in Sect. 14. 31 the 
subsequent developments are devoted to unitary weights 
W = Jjvf! they can easily be extended to include non 
unitary ones. Appendix iBl gives useful properties of these 
matrices. The reconstruction of x from y, i.e. the inversion 
of is a Fourier synthesis problem. 

In formulation © , the data y are in the (u, u)-plane 
while the map x is in the image plane. Two other state- 
ments are usually given: one regarding the (u, u)-plane 
only and the other the image plane exclusively. 

In terms of an usual approximation, after calibration, re- 
gridding, . . . Moreover, for the sake of readability, equations 
are given in ID and computation results are presented in 2D. 
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1. In the Fourier domain, © becomes a simple truncation 
by an invertible change of variable x = Fx: 

y = Tx + b. (2) 

Its inversion becomes a problem of extrapolating / in- 
terpolating "missing" Fourier coefficients. 

2. Furthermore, denoting y = T l y the zero-padded data, 
and y = F^y the dirty map, © becomes a convolution 

°y = Hx + b (3) 

where H = F^^TF is a (circulant) convolution 
matrix and b = T'6. (Superscripts "t" and "f" re- 
spectively denotes matrix transpose and conjugate- 
transpose). The instrument response (the dirty beam) 
is read in any one line of H, up to a circular shift. The 
inversion becomes a deconvolution problem. 

Remark 2. It should be noted, however, that the cor- 
relations of b and b differ from one another, and that 
in this sense, the two problems are not equivalent. 

Whichever formulation is envisaged, the tackled prob- 
lem is a rank-deficient linear inverse problem with additive 
noise. Indeed, the number of observed Fourier coefficient 
is far less than the number of pixels (M <C N) and the 
operators TF for Q, T for © or H for © have N — M 
singular values equal to 0, and M singular values equal 
to 1. Consequently, the least-squares criterion 

J™(x)=\\y-TFxf (4) 

possesses an infinite number of minimizers. The dirty map 
is one such solution since it cancels out J LS , and the other 
ones are obtained by adding maps with frequency compo- 
nents outside the (u, w)-plane coverage only. 

3. Regularization 

So, the selection of a unique solution requires a priori in- 
formation on the searched maps to be taken into account. 
In orde r to achieve this, we resort to regularization t ech- 
niques l|ldierll2001aHDemomentlll989HTarantolalll987l) . al- 
lowing diverse types of information to be considered, in 
order to exclude or avoid non desirable solutions. 

3.1. Criterion, penalization and constraints 

• Positivity and support. This information is naturally en- 
coded by hard constraints for the pixels. Let us note 
Ai, the collection of pixels on the map, S the collec- 
tion of pixels on a support, and S its complement in 
M. 

— (C s ) : support 

Vp € S , x p = . 

The proposed method takes into account the 
knowledge of a support (S is known and S ^ M) 
but remains also valid if S = M. 
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— (C p ) : positivity 

Vp e M , x p > o . 

This information is taken to be valid in the fol- 
lowing sections of this paper and all reconstructed 
maps will be positive. 

— (Ct) : template 

Vp G M , t~ <x p <t+ . 

It is also possible account for a known template but 
it is not numerically investigated in the paper. 
• Correlation structure. Here, we are concerned with the 
a priori correlation (ES) or non-correlation (PS) of the 
searched map. In the image plane, this information is 
naturally coded by penalization terms R(x), as a sum 
of potential functions </> which addresses the pixels. 

— (P c ) : the smooth map (ES) is favored by the in- 
troduction of interaction terms between pixels 

R c{x) = ^2<l>c [x q ,x p ] , (5) 

where p ~ q symbolizes neighbor pixels. 

— (P s ) : on the other hand, separable terms favor PS 

R s (x) =^0 s [x p ] . (6) 

These terms independently shrink the pixels to zero 
and therefore favor quasi-null maps. 

— (P m ) : in the following section, we will also need to 
penalize the average level of the maps 

R m {x) = <j> m [y^arj (7) 

so as to specifically compensate for the absence of 
Fourier coefficient at null-frequency. 

— (Pd) : it is also possible account for a known default 
map x through a specific penalization term such as 

R d (x) = ^2<f>d [x p ,x p ] (8) 
but this is not numerically investigated here. 

A criterion J is then introduced as a combination of 
some penalization terms ©-© and the data based one (@J 
according to the objective: PS component (Sect. l3~2"|l . ES 
component (Sect. li.-'A) and both of them simultaneously 
(Sect. GO}. In every case, the solution x is defined as the 
minimizer of J under constraints C p and C s : 

{min J(x) 
gt ix p = for peS (9) 
' 1 x p > for p e M 

that is to say, as the solution of problem (V) . One property 
then becomes crucial to the construction of J: 

— (Pi) : J is strictly convex and differentiable. 

Indeed, under this hypothesis, 



1. the problem (V) possesses a unique solution x, which 
allows the proper definition of the estimated map; 

2. the solution in question is continuous with respect to 
the data and to the tuning parameter values; 

3. a broad class of optimization algorithms is available. 

As J LS is itself (large sense) convex and differentiable, the 
property (Pi) can be assured if the potential functions are 
themselves convex and differentiable. Therefore, we resort 
to this kind of potential. 

Remark 3. Non-convex potentials hav e been introduced 
in image reconstruction in the 1980s l|Oeman fc Gemar] 
Il984t iBlake fc Zissermanlll987j) . As they are richer, they 
allow a sharper description of the searched images. For ex- 
ample, they can integrate binary variables, allowing con- 
tour detection to be carried out, at the same time as image 
reconstruction. As a counterpart, the involved criteria can 
possess numerous local minima. The computational cost 
for optimization then increases drastically, and sometimes 
without guarantee against local minima. 

In Sect.El several optimization schemes have been in- 
vestigated within the recommended convex framework. 
Various iterative algorithms solving (V) are concerned, all 
of them converging to the unique solution x whatever the 
initialization. The only question at stake is computation 
time. An other property of J is therefore crucial. 

— (P2) : J is quadratic and circular-symmetric. 

This property allows fast optimization algorithms to be 
put into practice taking advantage of the FFT algo- 
rithm: fast criterion calculations, explicit intermediate so- 
lutions,. . . Since J LS is itself quadratic and circulant, (P 2 ) 
is satisfied if the regularization terms are circulant and the 
potential functions <f> are Quadratic (Q) or Linear (L). 

Remark 4- Mixed convex potentials, generally quadratic 
about the origin and linear ab ove a certain threshold , 
are used in image processing l(Bouman fc Sauerl Il993h 
and especially in astronomical imaging ijMugnier et all 
in order to preserve possible edges. Fr om the opti- 
mization strategy stand point, recent works l)ldierll2001bl: 
lAllain et all2004[) allow to reduce the convex optimization 
problem to a partially quadratic one. This would make 
possible the development of an FFT and Lagrangian based 
algorithm for our PS+ES problem. We regard these forms 
as perspectives and we will see that forms Q and L are 
sufficiently rich and adapted to the envisaged contexts. 

3.2. Point Sources - separable linear penalty 



This section is devoted to PS: the proposed penalization 
term is of type © where <j) s is a potential of R + or W + 
onto R, to be specified. 



Classical MEM 1 


Nitvananda k. Naravanl 




1982} 


Naravan & Nitvananda 


1984: 


komesaroff et al 


1981: 


Gull & Daniell 1978: Gull & Skilling 1984: Bhandari 


1978: Le Bcsnerais et all Il999) 


come into play, when, 



for example, 4> s [x] = —log a;, <fi s [x] = xlogx or 



J.-F. Giovannelli and A. Coulais: Positive deconvolution for ES + SP. 



5 



4> s [x] = —x + x + x \ogx/x where x is a default 
map llO'SullivarJ 119951 ISnvder et al] Il992ft . They have 
been widely used in the domain and in ima ge recon- 
struction ijMohammad-Diaiari fc Demomentl Il988j) . 
They have the advantage of ensuring the property 
(Pi) on R^_, so the problem is properly regularized 
and (V) possesses a unique solution. They also enjoy 
the advantage of ensuring (strict) positivity, thanks 
to the presen ce of an infinite derivative a t the origin 
#(0+) = -oo frfamvan fc NitvanandallTflsi. 

However, these functions prohibit null-pixels and this 
can be seen as a flaw when the searched maps are largely 
made up of null-pixels. On the other hand, null-pixels are 
favored by the intr oduction of a potential </> s which pos- 
sesses at its origin llSoussenl lioOO') 

— a minimum value, and 

— a strictly positive derivative. 

Without loss of generality we set: S (O) = and </>s(0 + ) = 
1, while two possibilities allow property (P2) to be re- 
spected: the form L and the more general form Q. 
L : 4> s (x) = x 

Q : 4> s (x) — ax 2 + x 

The penalization is then written as: 



Rix) = A s x p + e s x 



(10) 



The strict convexity property (Pi) imposes e s > 0: the L 
term ensures a positive derivative at the origin, while 
the Q term ensures strict convexity. 

Remark 5. In order to favor high amplitude peaks, a least 
penalization function is desirable, i.e. e s = 0. In this case, 
it is possible that J remains unimodal or strictly convex, 
although we have no proof of this. This property could 
depend on the value of A s , on the knowledge and form of 
the support, on the (ii, i>) -plane coverage or on the data 
in each particular case. 



and even. In order to ensure property (P2), we are led to 
choose 4> c in class Q and to reject class L: cf) c (x) = x 2 and 

JV 

R(x) = A c ^2 \ x p+i ~ x pf 

with the hypothesis xq — xn in order to ensure circularity. 

We are here dealing with ear ly regularization tech- 
niques, that appeared in the 1960s l|Phillipsll962HTwomevl 
1963; lTikhonovlll963ft and were developed in the mid- 
1970s in works bv iTikho nov fc Arseninl ljl977|) in a con- 
tinuous context and bv lHuruTi|l977ft in a discrete context. 
They are also related to the well-known Wiener filter. 

In this form, the strict convexity condition (Pi) is not 
respected. Indeed, J LS is not sensitive to constant maps 
(since null-frequency is not observed) and neither is the 
regularization term (since it is only a function of the dif- 
ference between pixels). Several options are available for 
dealing with this indetermination. 

1. Support constraint C s : as soon as the support con- 
straint is valid, if at least one of the pixels is zero 
(S M), J is strictly convex on R 5 . 

2. In the absence of support information, it is sufficient 
to penalize the mean of the map by a term such: 



R, 



Intuitively, it reduces the mean of the map towards 
and is counterbalanced by the positivity constraint. 
3. It is also possible to penalize the quadratic norm of the 
map by a term such as that introduced in Sect. l3~2"l 

The penalization thus reads 

R(x) = A c ^ [ X P+l ~ x p} 2 + £ m \^2 x p\ • (11) 

Under this form, properties (Pi) and (P2) are satisfied if 
(e m > , A c > 0) in the case S = M and (e m > , A c > 0) 
in the case S ^ M.. 



3.3. Extended Sources - correlated quadratic penalty 

This section is devoted to ES: the penalty term of type © 
intr oduces interactio ns between neighboring pixels. 

lO'Sullivanl l)l995(l proposes the use of an I-divergence: 
4> c [x,x'\ = —x + x' + xlogx/x'or an Itakura-Saito dis- 
tance: cj> c [x, x'] = — logx/x' — 1 + x/x' in the symmetrized 
version. As in the case of Sect. 13.21 these allow property 
(Pi) and positivity to be ensured. However, they prohibit 
null-pixels and do not ensure property (P2)- 

We resort to classical terms of image processing based 
on finite differences between neighboring pixels. In the 
simplest case, first order differences yield 

4> c [x, x'] = (f) c [x- x'] 

where </> c is a potential of 1R onto 1R to be specified. In order 
to effectively favor smooth and correlated maps, and due 
to reasons of symmetry, <f> c is chosen to be minimal in 



3.4. Mixed model 

The present paper is devoted to maps composed of 
both type s of component sim ultan eously: ES and PS . 
Following l|Magain et al.lll998|) and l|Pirzkal et alJl2000ft . 
we introduce two maps x e and x p which describe each 
component respectively. The direct model Q becomes: 

y = TF(x c +x p )+b, (12) 

and the least-squares term 

J^l(x e ,x p ) - \\y-TF(x c + x p )\\ 2 , 

where the subscript "Mix" stand for Mixed map. This 
form raises new indeterminates as it now concerns the esti- 
mation of 2A^ variables, still from a single set of M Fourier 
coefficients. However, it allows the explicit introduction of 
characteristic information about each map through two 
adapted regularization terms. 
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1. A separable term for x p , identical to that in Sect. 13.21 

Rs(x p ) =J2x p (p), 

minimum at and with a strictly positive derivative. 

2. An interaction term between neighboring pixels of the 
map £c e , identical to that in Sect. EP1 

R c (x e ) = ^ i X e(p + !) - x e(p)] 2 ■ 

So as to ensure property (Pi), the same terms as in 
Sect. 13.21 and Sect. 13.31 are added, and the regularized cri- 
terion takes the form: 

■^Mjif (*"ej x p) = Jmi X ( x ct x p) (13) 
+ A c ^2 l x c(P + 1) - x c{p)f + Em Xe (P) ' 

where superscript "Reg" stands for Regularized. Regula- 
rization parameters (hyperparameters) A c and A s tune the 
smooth and spiky character of maps x e and x p . 

In this form, properties (Pi) — (P2) are satisfied if 
(A s > 0, A c > 0, e s > 0) and e m > when S ^ M or 
e m > when S = A4. The couple of maps (xl,x^) is 
properly defined as the solution of problem (V) and the 
next section (Sect.QJ gives the first practical results (sim- 
ulated and real data processing). Section |S] is devoted to 
a fast optimization algorithm. 

4. Computation results 

4.1. Nangay radioheliograph 

Radio emission of the Sun at meter wavelength is known 
since World War II. The Nancay radioheliograph (NRH) 
is a radio-interferometer dedicated to imaging the solar 
corona and it monitors the radio burst in solar atmosphere 
at such wavelengths with high temporal rate, adequate 
spatial resolution and high dynamic. 
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Fig. 2. Left figure shows instantaneous (u, u)-plane cov- 
erage (EW array is along vertical direction and NS array 
is along horizontal direction). Right figure gives the dirty 
beam, defined as the 2D Fourier transform of the (it, v)- 
plane coverage with a unitary weight for each visibility. 

At such frequencies, mainly two kinds of structures 
are observed in the corona: (1) larger structures (ES) 



and (2) smaller structures (PS). The quiet Sun is 
the largest structure, larger than the Sun size in the 
visible and slowly varying on long term scale (years) 
llLantos fc Alissa,ndrakisl \l99(j) . Medium size structures 
(1-m) are the radio co unterpart of coronal holes and mag- 
netic loops (plateau) ijAlissandrakis fc Lantodfl996|) . and 
are also observed simultaneously in soft X rays. The 
time scale for such structures is days to weeks. They 
are clearly correlated to persistent structures observed in 
other wavelength (optical and X-rays) and rotate on the 
radio maps quasi simultaneously with their optical and X- 
rays counterparts. The small structures (2) with very high 
brightness, can often reach several tens of Millions Kelvin 
l(Kerdraon fc MercierlllQSsTk they usually have a small life 
time (few seconds) and are associated to energetic events 
in the magnetic loops in the Sun's atmosphere. Correlation 
with structures observed in other wavelengths is more dif- 
ficult. 
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Fig. 3. Dirty maps typically encountered with NRH: sim- 
ulated data (top) and real data (bottom). Contour levels 
are — 10~ 2 to 5 10~ 2 , step 2.5 10 -4 (they are used for all 
the shown maps). 

The NRH is composed by two arrays: one along Est- 
West (EW) direction with 23 antennas, the other along 
North-South (NS) with 19 antennas. The NRH is op- 
erating in the range 150-450 MHz at a time sampling 
rate of 1/10 s., about eight hours a day, with favorable 
signal to noise ratio. Since the re furbishing of the in- 
strument l(Kerdraon fc Delouisll997^ cross-correlation be- 
tween most of the antennas in both arrays are available. 
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As a consequence: (i) 569 non redundant instantaneous 
visibilities are now available 2 (with unitary weights) and 
moreover, (ii) the instantaneous coverage of the (u, v)- 
plane (shown in Fig. |2J) becomes much more uniform. 
Nevertheless, due to the structure of the arrays, the cov- 
erage is not uniform. The central part of the (u, u)-plane 
essentially consists of two rectangular domains: the central 
one is a 16 x 16 square and the larger one is a 32 x 46 rect- 
angle. With this configuration 2D instantaneous imaging 
(without Earth rotation aperture synthesis) becomes pos- 
sible despite strong sidelobes in the dirty beam (see also 
Fig. |2l • As far as the dirty beam is concerned, the max- 
imum value is normalized to 1 and located in the middle 
of the map at (64, 64) . A secondary important lobe partly 
around (1, 64) and (128, 64) referred to as the aliasing lobe 
has amplitude 0.70 and characterizes important aliased 
response. The first negative lobe is —0.10 around the cen- 
tral lobe and —0.23 around the aliasing lobe. Moreover, 
the first positive lobe is 0.14 near the central lobe and 
0.12 near the aliasing lobe. In addition, the FWHM is 4.5 
(resp. 4) pixels for the central (resp. aliasing) lobe. 

At processed frequency (236 MHz), the field of view 3 
(FOV) related to the shortest baseline (55 m in NS, 50 m 
in EW) is ~ 1 20 ' and the size of the quiet sun is ~ 
40', i.e. ~ 1/2 FOV. Since at observed frequencies (150- 
450 MHz) the FWHM of the smallest antenna primary 
beam (few antennas are 15 m diameter) is much wider 
than the FOV, a unitary primary beam is appropriate. 
Moreover, the Shannon criterion is respected if the pixel 
number is - 60 for a FOV of 1 °. 

With the given characteristics, ES/PS separation must 
be achieved and reconstruction errors must be as small 
as possible for both maps, in order to strongly constrain 
physical models and to monitor position, amplitude and 
separation of bursts. But imaging the encountered con- 
text mixing PS and ES remains difficult and standard 
methods such as CLEAN and MEM (even in a multireso- 
lution approach) are usually inefficient due to the large 
background a nd the intrica te mixing of real structures 
and sidelobes l|Coulaisll997(l . One possible outcome of the 
present work is to provide to the solar radio community 
accurate maps from NRH in order to achieve more de- 
tailed scientific studies. The following computation study 
(simulated and real data) is a typical case encountered 
with NRH and provides a first element in this sense. 

4.2. Simulation results 

Simulated data 

The true ES map x* (Fig.^-a), is ranging in amplitude 
(arbitrary units) from to 5.5 10 -3 . The true Sun lies in 

2 Thanks to Hermitian symmetry 1138 Fourier coefficients 
are available. The computed (u, v)-p\a,ne and map are 128x128. 

3 For a declination 23 ° and null hour angle (the Sun at noon 
in summer), the FOV is 87' in EW and 86' in NS, and the 
resolution is 3.27' in EW and 2.17' in NS, since main EW arm 
is 1600 m with step 50 m and NS arm is 2640 m with step 55 m. 



a disk centered in the middle of the image, i.e. (64,64) 
with a 64 pixels diameter. The outer part of the disk is 
zero and the mean of this component is 5.59 10~ 4 . The 
true PS map x* (Fig.0}b) consists of two peaks: the first 
one is located at (60,61) with amplitude 5.0 10~ 2 , and 
the second one overlaps pixels (57, 61) and (57, 62) with 
respective amplitudes 5.0 10~ 2 and 4.5 10 -2 . Data (in the 
(u, u)-plane) are simulated using the direct model Ea. lfT^I . 
i.e. FFT and truncature, and corrupted by a white, zero- 
mean complex Gaussian noise with variance 2 10 -7 . (This 
noise variance has been chosen in order to mimic real 
data). The dirty map is shown in Fig. EI If is clearly dom- 
inated by the PS, and the whole map is corrupted by side 
lobes. Moreover, the two close peaks at location (60,61) 
and (57,61) — (57,62) are not resolved. 

Reconstruction parameters 

The supports have been deduced from the dirty map. It 
is a disk centered at (64, 64) with a 70 pixels diameter for 
the ES map. Regarding the PS map, the support consists 
of one disk centered in (58, 61) with a 10 pixels diameter. 

In practice, two hyperparameters have to be tuned: A c 
and A s (e s is practically set to 10~ 10 ). A c must be set in 
the order of magnitude of eigenvalues of the Hessian of 
the criterion and is set to A c = 2. A s has been empirically 
selected after several trials in order to visually achieve 
separation of PS and ES: it has been set to A s = 10~ 3 . 

Reconstruction results 

Fig.0]shows the reconstructed maps. A simple qualita- 
tive comparison with the references x c and x p shows that 
the two components x~ a and x^ are efficiently separated 
and accurately reconstructed. 

The two peaks of x^ shown in Fig. 0J-d are precisely 
located at (60,61) and (57,61) — (57,62) (overlapping). 
The estimated amplitudes are 0.051, 0.048 and 0.043 re- 
spectively, i.e. an error of less than 5%. Moreover, the 
two close peaks are separated whereas they are not in the 
dirty map. This illustrates the resolution capability of the 
proposed method resulting from both data and accounted 
information (positivity, support, and PS+ES hypothesis). 
It is also noticeable that the respective part of flux in 
overlapped pixels (57, 61) — (57, 62) is correctly restored. 

Fig. 2Jc gives the estimated ES map xl. Compared to 
the true one of Fig. QJa the main structures are accurately 
estimated. The contour lines of Fig. 0J-C are very similar 
to the one of Fig.^Ja and the relative reconstruction error 
is less than 2%. Moreover, the mean of the estimated ES 
map x^ is 5.57 10 -4 while the true mean is 5.59 10 -4 : the 
total flux is correctly estimated. The maximum value is 
5.4 10 -3 in xl whereas it is 5.5 10 -3 in x c : the dynamic 
is also correctly retrieved. Nevertheless, a slight distor- 
tion located around pixel (65, 60) can be observed in the 
proposed ES map. It probably results from an imperfect 
separation of the two components: a slight trace of the 
dirty beam remains in the estimated ES map. Moreover, 
the sharp edges of the true Sun are slightly smoothed due 
to the lack of high frequencies in the available Fourier 
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a: True object x% 



20 40 60 80 100 120 

b: True object x* v 




20 40 60 80 100 120 

c: Estimated object x2 



20 40 60 80 100 120 

d: Estimated object x~p 



Fig. 4. Simulation results (see Sect. I4.2|> . Contour levels are the same than in Fig. and for all the maps: the true 
ES x* (a) and the estimated one x^ (c) as well as for the true PS x* (b) and the estimated one x^ (d). 




a: Estimated object xZ 



b: Estimated object x p 



Fig. 5. NRH data processing from typical scientific observation at 236 MHz (see Sect. I4.3|) . Contour levels are the 
same than in Fig. [3] and QJ The two components x~l (a) and x~p (b) are clearly separated and deconvolution of both 
component is clearly achieved. Both maps are positive and the prescribed supports are respected. 



coefficients incompletely enhanced by accounted prior in- 
formation (see Remark and 0J) . 

4.3. Real data computations 

This section is devoted to real data processing based on a 
data set from the NRH 4 . The coverage is identical to the 
one of simulated data of the previous Section. 

4 The eleventh of June, 2004, 13h00, at 236 MHz. 



The dirty beam is shown in Fig. [21 and the dirty map is 
shown in Fig. [3J Both dirty beam and dirty map are typ- 
ically encountered with NRH and are similar to the one 
simulated in the previous section. As expected, resolution 
is limited and the quality of the map is entirely contin- 
gented upon sidelobes around the brightest point sources 
(radio burst). Imaging such a complex context mixing PS 
and ES suffers from intricate mixing of real structures and 
sidelobes due the brightest ones. 
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The same supports have been used to compute the real 
data and the simulated ones. It is a disk centered in the 
middle of the map with a 70 pixels diameter for the ES 
map and a disk centered in (58, 61) with 10 pixels diameter 
for the PS map. The same value of the parameters A c = 2 
and A s = 10 -3 have been used to compute the real data 
and the simulated one (e s remains set to 10 -10 ). 

Estimated maps are shown in Fig. [SJ-a (ES compo- 
nent) and Fig. ISJ-b (PS component). The two components 
xl (Fig. |SJ-a) and x~p (Fig. [5Jb) are clearly separated and 
both deconvolution is clearly achieved. Both maps are pos- 
itive and the prescribed supports are respected. Moreover, 
the map presents a similar structure to the usual one 
of the Sun at meter wavelengths without strong point 
sources ()Coulaislll997t Ihantos fc Alissandrakislll996|) . 

5. Numerical optimization stage 

The estimated maps are defined as the unique solution of 
the problem (V) given by (jSJ) which involves the quadratic 
criterion J given by i|13f) . Up to an additive constant: 



penalty term cXp/2. The entire term write: 



J(x) 



q t x 



where ] collects the two maps (appendix O 

gives Q and q). Thus, (V) is a convex quadratic program: 



min — x t Q x + q* x 



(V) I 



s.t. 



x p = for p G S 
x p >0 for p G M 



(14) 



widely investigated in the optimization literature. The 
main difficulty is twofold. On the one hand, the non- 
separability of J together with positivity constraint pre- 
vents from explicit optimization. On the other hand, the 
number of variables is very large. We have investigated 
most of the proposed methods in the excellent reference 
book (iNocedal fc Wriehtll2000h : 

— Constrained gradient 

— Gradient projection 

— Barrier and interior point 

— Relaxation (coordinate-by-coordinate) 

— Augmented Lagrangian (method of multipliers) 

and have selected the latter as the faster. It is based 
upon successive optimizations of a Lagrangian function C 
founded on Lagrange multipliers £, slack variables s and 
quadratic penalty. It is computationally based on FFT and 
threshold, so it is, in addition, very simple to implement. 



£p x p + r, c x p ■ 



(15) 



pes 



pes 



The inequality constraint x v > (p £ S) is converted 
into the equality one s p — x p = using the slack variable 
s p > 0. Lagrange and penalty terms then write: 



£p(x p S p ) + C ^""^(Xp Sp) 



(16) 



pes 



pes 



In order to simultaneously process both equality l|15fl 
and inequality l|16(l constraints, we introduce extra slack 
variables s p — for p£5. The Lagrangian then writes: 

C(x, s,£) = J(x) - t(x - s) + - c(x - sf(x - s) 
where s and £ collect slack variables s p and multipliers £ p . 



5.2. Algorithm 

The algorithm then iterates three steps: 

(T) unconstrained minimization of C w.r.t. x, 

(2) minimization of C w.r.t. s, s.t. s p > 0, 

(3) update t and c. 

The efficiency of the proposed algorithm relies on both 
slack variables and property (P2)- Roughly speaking, pos- 
itivity is transfered on slack variables, so, the non separa- 
ble constrained problem (V) is split in two subproblems: a 
non-separable but unconstrained one computable by FFT 
(step ©) and a constrained but separable one (step (2)). 

Step (1) proceeds by fixing t and s to the current value 
and then computes the unconstrained minimizer x of C. 
It is an unconstrained convex quadratic problem, so its 
solution is explicit: 

x = -{Q + c/jv) -1 (q + [£ + cs}) , 

and computable by means of FFT, thanks to circularity. 

Step @ updates the slack variables s p for p 6 S (by 
construction, s p = for p £ S) as the minimizer 7s p of £, 
subject to s p > 0. 



max (0, cx p — £ p )/c 




for p G S 
for p G S 



This step is constrained but separable: the constrained 
minimizer is the unconstrained one if positive and if not. 
Step @ consists in updating the Lagrange multiplier £: 



max (0, £ p — cx p ) 



for p G S 
for p G S 



5.1. Lagrangian Function 

The equality constraint x p = (p G S) is introduced by 
means of a usual Lagrangian term —£ p x p together with a 



This step can also include an update of c (e.g. c = 1.1c). 
Practically, c is not updated (see next subsection). 

Steps (1) to (3) are iterated until stopping condition is 
met, e.g. relative variation smaller that 0.1%. 
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Remark 6. Constrained variables x p — for p € S can 
also be eliminated. This is a relevant strategy when using 
gradient based or relaxation methods. It does not prevent 
from computing J and its gradient by means of FFT. On 
the contrary, such a strategy is not relevant here: it would 
break circularity and prevent from using FFT in step (T). 

5.3. Practical case and computations time 

This section specializes the algorithm in the case of con- 
stant coefficient c. In this case, step ©-© reduces to: 



for peS 
for p £ S 



Moreover, Q + cIn can be inverted once for all and the al- 
gorithm then requires 4 FFT per iteration. The algorithm 
has been used in the previous computations with constant 
coefficient c = 10 -3 . Convergence is achieved after about 
1000 iterations and it takes half a minute 5 . 



6. Conclusions 

The problem of incomplete Fourier inversion is addressed 
as it arises in map reconstruction (deconvolution, spectral 
interpolation/extrapolation, Fourier synthesis). The pro- 
posed solution is dedicated to specific situations in which 
the imaged object involves two components: (i) an ex- 
tended component together with (ii) a set of point sources. 
For this ca ses, new development s are given based on exist - 
ing work of lMagain etaD jl99Sh and lPirzkal et all i|2000t . 

The main part of the paper deals with inversion in the 
regularization framework. It essentially departs from usual 
strategies by the way it accounts for (1) noise and inde- 
terminacies, (2) smoothness/sharpness prior and (3) posi- 
tivity and support, in a unique coherent setting. The pre- 
sented development can also include known template and 
default map. Thus, a new regularized criterion is intro- 
duced and estimated maps are properly defined as its 
unique minimizer. The criterion is iteratively minimized 
by means of an efficient algorithm essentially based on 
Lagrange multiplier which practically requires FFT and 
threshold only. The minimizer is shown to be both prac- 
tically reachable and accurate. A first evaluation of the 
proposed method has been carried out using simulated 
and real data sets. We demonstrate ability to separate the 
two components, high resolution capability and high qual- 
ity of each map. To our knowledge, such a development is 
an original contribution to the field of deconvolution. 

Nevertheless, a further evaluation of the proposed 
method is desirable. Future work will include systematic 
evaluation of the capability of the proposed method as 
a function of (u, w)-plane coverage, PS amplitudes versus 
ES ones, PS position (especially in a subpixelic sense) and 



5 Algorithm has been implemented with the computing envi- 
ronments Matlab and IDL on a PC, with a 2 GHz AMD- Athlon 
CPU, and 512 MB of RAM. Both codes are ~ 50 lines long. 



noise level. Such an assessment concerns both simulated 
and real data. Moreover evaluation of the potentiality of 
the method on large maps (e.g. VLA images), high dy- 
namic imaging (e.g. WSRT images) and imaging using 
millimeter interferometers (e.g. IRAM PdBi and ALMA) 
or optical instruments will be considered. 

A part of future work in the field of SP+ES imaging, 
will include convex non quadratic penalization of ES (see 
Remark 0J). Another part of future work will particularize 
the proposed method in order to produce maps of ES only 
and maps of PS only. 

A Bayesian interpretation of the proposed method in- 
volve truncated Gauss-Markov models (ES component) 
and exponential white noise (PS component) and formally 
provides likelihood tools in order to achieve automatic 
tuning of the hyperparameters. This is a more delicate 
aspect but it will be addressed in future works. 
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Appendix A: Notations 

In this paper, Ip denotes the P x P identity matrix and 
(resp. Af*) denotes the complex conjugate transpose 
(resp. transpose) of a given matrix M. 

Let us note D the (circulant) first order difference ma- 
trix and Ad = FDF^ the diagonalized matrix. Let us 
also note 1 the ones column vector with TV components 

o 

and 1 — Fl its FFT (non- null at null frequency only). 

We introduce now two matrices Ae and Ap useful to 
compute ES and PS respectively: 



oto 

1 1 



where Ac = T l T. The three sub-matrices Ae, Ap and 
Ac arc diagonal matrices. 

Appendix B: Conventions and properties 

This Appendix gives several properties of F and T intro- 
duced in Sect. 

— F* F = FF^ = In ■ orthonormality of the normalized 
FFT. 

— T is a truncation operator, N x M, (eliminates coeffi- 
cient outside the coverage). 

— T* is a zero-padding operator, M x N, (adds null co- 
efficient outside the coverage). 

— Ac = T*T is a projection matrix, N x N, (nullifies 
coefficients outside the coverage). 

-TT t =I M . 
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Table C.l. Functions, gradients and Hessian of encountered criteria (given as a function of map and their FFT). 



Appendix C: Gradient and Hessian calculi 



The update reads : 



This Appendix is devoted to the vector q and the matrix 
Q involved in the minimized criterion. 

q is a 2N components column vector based on the 
gradient of criterion J, at x = 0. The first part is the dirty 
map and the second one is the dirty map minus a constant 
map equal to A s /2. In the Fourier domain, q reads: 



o _ dJ_ 
ox 



x=0 



dJ 
dx e 

dJ 
dx D 



= -2 



y 

y A s i/2 



x=0 



Q is a 2N x 2N matrix based on the Hessian of J. The 
two anti-diagonal elements arc the Hessian of the LS term 
and rely on the dirty beam only. The diagonal elements 
are the Hessian of J w.r.t. each map x p and x c . In the 
Fourier domain, Q reads: 



° d 2 J 
ox 



d 2 j 

r\ O 2 

Ox c 



d 2 j 

r\ O <~| O 

OXp ox c 



d 2 J 



OX c OXr 



d J 



dxr 



A E 
A c 



A c 
A P 



Appendix D: Object updates 

The present subsection gives details about the step of 
the proposed algorithm (Sect.0) : the unconstrained min- 
imization of C w.r.t. x, i.e. the update of x e and x p . Let 
us introduce the two vectors 

(z c = y+{i c + c°s c )/2 

\z p = y+{£ p + cs p )/2-X s l 

based on observed data y and FFT of slack variables and 

o ° 

Lagrange multipliers s = Fs and £ = F£ (for each map 
ES and PS). Let us also introduce two diagonal matrices 



M E = A E 
M P = A P 



cI N /2 
cI N /2 



x c = (M E M P - A c ) (M P z c - A c z p ) 
x p = (Me Af p - A 2 ,)' 1 (M E z p - A c z c ) 



easily implemented since Me Mp — A c is diagonal. 
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